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Abstract 

Cn ' We present an extension of our GPGCD method, an iterative metliod for cal- 

culating approximate greatest common divisor (GCD) of univariate polynomials, to 
polynomials with the complex coefficients. For a given pair of polynomials and a de- 
gree, our algorithm finds a pair of polynomials which has a GCD of the given degree 
and whose coefficients are perturbed from those in the original inputs, making the 
perturbations as small as possible, along with the GCD. In our GPGCD method, the 
problem of approximate GCD is transfered to a constrained minimization problem, 
r ) , then solved with a so-called modified Newton method, which is a generalization of the 

^^ ' gradient-projection method, by searching the solution iteratively. While our original 

^^ ' method is designed for polynomials with the real coefficients, we extend it to accept 

(~| . polynomials with the complex coefficients in this paper. 

-4— » 

1 Introduction 

For algebraic computations on polynomials and matrices, approximate algebraic algorithms 
^ ' are attracting broad range of attentions recently. These algorithms take inputs with some 

^T , "noise" such as polynomials with floating-point number coefhcients with rounding errors, or 

^^ ' more practical errors such as measurement errors, then, with minimal changes on the inputs, 

seek a meaningful answer that reflect desired property of the input, such as a common factor 
of a given degree. By this characteristic, approximate algebraic algorithms are expected 
to be applicable to more wide range of problems, especially those to which exact algebraic 
algorithms were not applicable. 

As an approximate algebraic algorithm, we consider calculating the approximate greatest 
common divisor (GCD) of univariate polynomials, such that, for a given pair of polynomials 
and a degree d, finding a pair of polynomials which has a GCD of degree d and whose coef- 
J^ ' ficients are perturbations from those in the original inputs, with making the perturbations 

H I as small as possible, along with the GCD. This problem has been extensively studied with 

various approaches including the Euclidean method on the polynomial remainder sequence 
(PRS) ([1], [15, m]), the singular value decomposition (SVD) of the Sylvester matrix 
([3] I E]): the QR factorization of the Sylvester matrix or its displacements ([1], [H], pU]). 
Pade approximation TT, optimization strategies ([2], |7], [5], [5], [H]). Furthermore, stable 
methods for ill-conditioned problems have been discussed (|3], [TU], [T5]). 

Among methods in the above, we focus our attention on optimization strategies. Already 
proposed algorithms utilize iterative methods including the Levenberg-Marquardt method 
[2], the Gauss- Newton method 19 and the structured total least norm (STLN) method 
(H) [H])- Among them, STLN-based methods have shown good performance calculating 
approximate GCD with sufficiently small perturbations efficiently. 
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In this paper, we discuss an extension of the GPGCD method, proposed by the present 
author [17] . an iterative method with transferring the original approximate GCD problem 
into a constrained optimization problem, then solving it by a so-called modified Newton 
method [IH] , which is a generalization of the gradient-projection method [12] . In the pre- 
vious paper 1171 . we have shown that our method calculates approximate GCD with per- 
turbations as small as those calculated by the STLN-based methods and with significantly 
better efficiency than theirs. While our original method accepts polynomials with the real 
coefficients as inputs and outputs in the previous paper, we extend it to handle polynomials 
with the complex coefficients in more generalized settings in this paper. 

The rest part of the paper is organized as follows. In Section [2j we transform the 
approximate GCD problem into a constrained minimization problem for the case with the 
complex coefficients. In Section [31 we show details for calculating the approximate GCD, 
with discussing issues in minimizations. In Section 21 we demonstrate performance of our 
algorithm with experiments. 

2 Formulation of the Approximate GCD Problem 

Let F{x) and G{x) be univariate polynomials of degree m and n, respectively, with the 
complex coefficients and m > n > 0. We permit F and G be relatively prime in general. 
For a given integer d satisfying n > d > 0, let us calculate a deformation of F{x) and G{x) 
in the form of 

F{x) = F{x) + AF{x) = H{x) ■ F{x), G{x) = G{x) + AG{x) = H{x) ■ G{x), 

where AF{x) and AG{x) are polynomials with the complex coefficients, whose degrees do 
not exceed those of F(a;) and G{x), respectively, H{x) is a polynomial of degree d, and F{x) 
and G{x) are pairwise relatively prime. In this situation, H{x) is an approximate GCD of 
F{x) and G{x). For a given d, we try to minimize ||/ii^(x)||2 -1- ||ZiG'(a;)||2 the norm of the 
deformations. 

In the case F{x) and G{x) have a GCD of degree d, then the theory of subresultant tells 
us that the ((i—l)-th subresultant of F and Gbecomes zero, namely we have Sd_i(i^, G) = 0, 
where Sfc(F, G) denotes the fc-th subresultant of F and G. Then, the (d— l)-th subresultant 
matrix Nd-i{F,G), where the fc-th subresultant matrix Nk{F,G) is a submatrix of the 
Sylvester matrix N{F, G) by taking the left n — k columns of coefficients of F and the left 
m — k columns of coefficients of G, has a kernel of dimension equal to 1. Thus, there exist 
polynomials A{x),B{x) G C[a;] satisfying 

AF + BG = 0, (1) 

with deg(A) < n—d and deg(i?) < m — d and A{x) and B{x) are relatively prime. Therefore, 
for the given F{x), G{x) and d, our problem is to find AF{x), AG{x), A{x) and B{x) 
satisfying Eq. ([Ij with making ||zlF||2 -I- ||^G||2 as small as possible. 
Let us assume that F{x) and G(x) are represented as 

F{x) = (/™4 + f^^2i)x"' + ■■■ + (/o,i + /o,2i) = Frc(x) + iFu^ix), 
G{x) = (.g„,i + gn,2i)x'' H h (.go,i + .90,2*) = Grc{x) + iGim{x), 

where fj^i, gj^i, fj^2, gj,2 are the real numbers and i is the imaginary unit, and F^dx) and 
GRe(a^) represent the real part of F{x) and G{x), respectively, while i^im(2;) and Gim(x) 



represent the imaginary part of F(x) and G{x), respectively. Furthermore, we represent 
F{x)^ G{x), A{x) and B{x) with the complex coefficients as 

Fix) = ifraA + fmai)x"' + ■■■ + (/o,l + fQ,2i)x° = Fr,c(x) + iPi^ix), 

G{x) = (g„4 + gn.2i)x" H h {gox° + go,2i)x" = GRo(a;) + iGim{x), 

A{x) = (a„_d,i + a„_d,2«)a;""'* H h (ao,i + aQ^2i)x° = v4Rc(a;) + iAim(a;), 

B{x) = (6™_d,i + fom-d,2i)a;™"'^ + • • ■ + (60,1 + 6o,2i)a;° = Bnc{x) + iBi^{x), 

respectively, where fj^i, fj^2, gj,i, 5j,2, aj,i, clj,2, &j,i, &j,2 are the real numbers, and, as in 
above, F^dx), Gyic{x), A^ie^x) and i?Re(a;) represent the real part of F(a;), G{x), A{x) and 
B{x), respectively, while -Fi„i(a;), Gimix), Aiai{x) and Biai{x) represent the imaginary part 
of -F(x), G{x), A{x) and B{x), respectively. 

For the objective function, ||Zii^||| + ||^G||2 becomes as 

m n 

E[(/^M - /^m)' + (/^-^ - h^f] + EtfeM ' 9,sf + {h2 - 9,.2n (3) 

For the constraint, Eq. ([Ij becomes as 

Nd-i{F, G) ■ * {an-dA + an-d,2i, • • • , ao,i + ao,2«, 6m-d,i + Ki-dah ■ ■ ■ , &o,i + ^0,2*) 

= 0. (4) 

By expressing the subresultant matrix and the column vector in Q separated into the real 
and the complex parts, respectively, we express ^ as 

{Ni+N2i)ivi+V2i) = 0, (5) 

Ni - Nd-i{Fjie{x),GR,{x)), N2 = Nd-i{Fi,,,{x),Gi^{x)), 

t t ^ ' 

Vl — {an-d,lT ■ ■ ,ao,l,bm-d,l, ■ ■ ■ ,bos), V2 = {an-d,2, ■■■ ,0-0,2, bm-d,2, ■■■ ,bo^2)- 

We can expand the left-hand-side of Eq. ^ as (A^i -I- N2i){vi + V2i) = {NiVi — N2V2) + 
i{NiV2 + N2V1), thus, Eq. ([5]) is equivalent to a system of equations: NiVi — A^2'U2 = 
0,NiV2 + N2V1 = 0, which is expressed as 

m -N2\(v,\ 

N2 N, J yv2j ''■ ^'> 

Furthermore, we add another constraint for the coefficient of A{x) and B{x) as 
\\A{x)\\l + \\B{x)\\l = K-d,! + ■■■+ «o,i) + (bl-d,i + ■■■ + bli) 

+ {al-da + ■■■ + «o,2) + {bL-da + ■■■ + ^0,2) -1 = 0, (8) 

which can be expressed together with ([7]) as 

-1\ /^i\ 

= 0, (9) 

where Eq. ^ has been put on the top of Eq. ([7]). Note that, in Eq. ([9]), we have total of 
2{m + n — d+l) + \ equations in the coefficients of polynomials in ([5]) as a constraint, with 
the j'-th row of which is expressed as qj — Q. 
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Now, we substitute the variables 

{/m,l' ■ • ■ j/o,lj5n,l: • • • 7 90,1, fm,2, ■ ■ ■ ,fo,2,gn,2, ■ ■ ■ ,ff0,2, 

fln-d,!, • • ■ 7 ao,l; b„i^d,l, ■ ■ ■ , bo^i,an-d,2, ■ ■ ■ , 0,0^2, bm-d,2, ■ • • , &0,2), (10) 

as a; = (xi, . . . , X4(m+n-d+2))i then Eq. ([3]) and ([9]) become as 

f{x) = {xi - /r„,i)^ H h (a;m+i - /o,i)^ + (2;m+2 - 9n.if H h (a;TO+„+2 - go,i)^ 

+ (a^m+n+S -~ fm,2) + ■ ' ■ + (x2m+n+3 ^ /o,2) 

+ {X2m+n+i - gn,2) + 1" (a;2(m+n+2) ^ 50,2) , (H) 

q{x) ^ \qiix), . . . ,q2(rn+n-d+l)+lix)) =0, (12) 

respectively. Therefore, the problem of finding an approximate GCD can be formulated as 
a constrained minimization problem of finding a minimizer of the objective function f{x) 
in Eq. (HI]), subject to q{x) = in Eq. 1^^. 



3 The Algorithm for Approximate GCD 

We calculate an approximate GCD by solving the constrained minimization problem (jlip . 
(IT^ with the gradient projection method by Rosen [T^] (whose initials become the name 
of our GPGCD method) or a modified Newton method by Tanabe [TB] (for review, see the 
author's previous paper [U]). Our preceding experiments [T71 Section 5.1] have shown that 
a modified Newton method was more efficient than the original gradient projection method 
while the both methods have shown almost the same convergence property, thus we adopt 
a modified Newton method in this paper. 

In applying a modified Newton method to the approximate GCD problem, we discuss 
issues in the construction of the algorithm in detail, such as 

• Representation of the Jacobian matrix Jg{x) and certifying that Jg{x) has full rank 
(Section EH]), 

• Setting the initial values f Section 13. 2p . 

• Regarding the minimization problem as the minimum distance problem (Section l3.3p . 

• Calculating the actual GCD and correcting the coefficients of F and G f Section 15^ . 
as follows. 



3.1 Representation and the rank of the Jacobian Matrix 

Pnx" -\ h poa;°, let Ck{P) be a 



For a polynomial P(x) G C[a;] represented as P{x) 
complex (n + fc, fc + 1) matrix defined as 
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For co-factors A{x) and B{x) in ^, define matrices Ai and A2 as 

Al = [C„,iAncix)) CniB^eix))], A2 = [CrrMlmix)) C„ (Bi„,(a:))] . (13) 

(Note that Ai and A2 are matrices of the real numbers oi m + n — d-\-l rows and m + n + 2 
columns.) Then, by the definition of the constraint (J12p . we have the Jacobian matrix Jg{x) 
(with the original notation of variables (jlO|) for x) as 

/ 2 • *t;i 2 • *V2^ 

Jgix) = Al -A2 iVl ~N2 

\A2 Al N2 Ni 

with Al and A2 as in p^ and A^i, iV2, fi and V2 as in ([5]), respectively, which can be easily 
constructed in every iteration. 

In executing iterations, we need to keep that Jg{x) has full rank: otherwise, we are 
unable to decide proper search direction. For this requirement, we have the following 
observations. 

Proposition 1. Let x* £ Vg be any feasible point satisfying Eq. (|12p . Then, if the cor- 
responding polynomials do not have a GCD whose degree exceeds d, then Jg{x*) has full 
rank. 

Proof. Let X* = {fm^l, ■ . . , fo,l,gn,l, ■ ■ ■ ,go,l, fm,2, ■ ■ ■ , fo,2,gn,2, ■ ■ ■ ,go,2,an-d,l, ■ ■ ■ ,ao,l, 

bm-d,i,- ■ ■ ,bo^i,a„-d,2,- ■ ■ ,ao^2,bm-d,2,- ■■,bo,2) with its polynomial representation 
expressed as in ^ (note that this assumption permits the polynomials F{x) and G{x) to 
be relatively prime in general). To verify our claim, we show that we have rank(Jg(a;*)) = 
2(m + n — d + 1) + 1. Let us express Jg{x*) — [Jh \ Jr), where Jl and Jr are column 
blocks expressed as 

/2 -'vi 2 • *t;2^ 
Jr = A^i -N2 
\ N2 Ni 

respectively. Then, we have the following lemma. 
Lemma 1. We have rank( Jl) = 2(m + n — d + 1). 




Proof. F or Ai = [Cm{A)i C„(i?)i], let Cm{A)i be the right m — d columns of Cm(A)i and 
Cn{B)i be the right n — d columns of C„(i?)i. Then, we see that the bottom m + n — 2d 
rows of the matrix C — [Cm{A)i C„(i3)i] is equal to the matrix consisting of the real part 
of the elements of N{A,B), the Sylvester matrix of A{x) and B{x). By the assumption, 
polynomials A{x) and B{x) are relatively prime, and there exist no nonzero elements in C 
except for the bottom m + n — 2d rows, thus we have rank(C) ^ m + n — 2d. 

By the structure of C and the lower triangular structure of Cm{A)i and C„(i3)i, we 
can take the left d + 1 columns of CmiA)i or C„(i?)i satisfying linear independence along 
with C, which implies that there exist a nonsingular square matrix T of order m + n + 2 
satisfying 

AiT = R, (14) 

where i? is a lower triangular matrix, thus we have rank(yli) = rank(i?) = m + ri — d + 1. 



Furthermore, by using T and R in (J14p . we have 
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foUowed by a suitable transformation on columns on the matrix in the right-hand-side of 
(lisp , we can make A2T to zero matrix, which implies that 



rank(JL) = rank R -A2T = 2 • rank(i?) = 2(m + n - d + 1). 









R 


-A2T 


hT 


R 



This proves the lemma. D 

Proof of Proposition [T] (continued). By the assumptions, we have at least one 
nonzero coordinate in the top row in Jr , while we have no nonzero coordinate in the top row 
in Jl, thus we have rank( Jg(a;)) = 2{m -t- n — d -|- 1) -I- 1, which proves the proposition. D 

Proposition [T] says that, so long as the search direction in the minimization problem 
satisfies that corresponding polynomials have a GCD of degree not exceeding d, then Jg{x) 
has full rank, thus we can safely calculate the next search direction for approximate GCD. 

3.2 Setting the Initial Values 

At the beginning of iterations, we give the initial value Xq by using the singular value decom- 
position (SVD) of iV = (^'^ ~^A in as iV = (7 S *F, [/ = (wi, . . . , M2(m+„-2d+2)), 

S = diag{ai,...,a2(m+n-2d+2)),V = (Vl, . . . ,t'2(m+n-2d+2)), with Uj e Jl2{m+n-d+l) ^ 

Vj £ jj^2(m+n-2d+2)^ ^^^ j^ _ (jia,g(cri , . . . , cr2(m+n-2d+2) ) dcuotcs the diagonal matrix 
whose the j-th diagonal element is aj. Note that U and V are orthogonal matrices. Then, 
by a property of the SVD O Theorem 3.3], the smallest singular value cr2{m+n-2d+2) 
gives the minimum distance of the image of the unit sphere s2(m-i-n-2 +2)-i^ given as 
g2(m+«-2d+2)-i ^ {a. £ ^2{m+n-2d+2) | ||^||_^ ^ ^|^ ^^ ^^ represented as 

with a2{m+n-2d+2)U2{m+n-2d+2) as its coordiuatcs. Thus, we have N_ ■ f2(m+n-2d+2) = 

0'2(m+n-2d-|-2)W2(m+™-2d+2)- For Vm+n-2d — {o-n-d, • • ■ , So, hn-d, ■ ■ ■ , &o), let 

A{x) = a„_dx"~'' + • ■ • + a()X° and B{x) = 5„j_d2;'""'* + • ■ • + b()X° . Then, A(x) _and 
B{x) give the least norm of AF + BG satisfying ||A||2 + ||i3|l2 = 1 by putting A{x) = A{x) 
and B{x) = B{x) in ©. 

Therefore, we admit the coefficients of F, G, A and B as the initial values of the 
iterations as 

Xo = (/m.l, • ■ • , /o,l,ffn,l, ■ • ■ ,30.1, /m,2, ■ • ■ ,/o,2,5n,2, • ■ • ,50,2, 

0.n-d,l, ■ ■ ■ , flO.lj bn-d,l, ■ ■ ■ , &0,1: 0'n-d,2, ■ ■ ■ , 00,2, &n-d,2, • ■ • , ^0,2)- 



3.3 Regarding the Minimization Problem as the Minimum Dis- 
tance (Least Squares) Problem 

Since we have the object function / as in dill) , we have 

V/(a;) = 2 • {xi — fm,l, ■ ■ ■ ,Xm+l — foS,Xm+2 — Qn,!, ■ ■ ■ ,Xm+n+2 ~ 50,1, 

Xm+n+Z ~ fm,2, ■ ■ ■ , X2m+n+3 ~ /o,2, X2m+n+4 ~ 9n,2, ■ ■ ■ , a^2(m+n+2) ~ ff0,2, 0, . . . , 0). 

However, we can regard our problem as finding a point x E Vg which has the minimum 
distance to the initial point Xq with respect to the (xi, . . . ,X2(„i+ri+2))"COordinates which 
correspond to the coefficients in F{x) and G{x). Therefore, as in the real case (see the 
authors previous paper [f 7'), we change the objective function as f{x) = :^f{x), then solve 
the minimization problem of f{x), subject to q{x) = 0. 

3.4 Calculating the Actual GCD and Correcting the Deformed 
Polynomials 

After successful end of the iterations, we obtain the coefficients of F{x), G{x), A{x) and 
B{x) satisfying ([1]) with A{x) and B(x) are relatively prime. Then, we need to compute 
the actual GCD H{x) of F{x) and G{x). Although H can be calculated as the quotient of 
F divided by J5 or G divided by A, naive polynomial division may cause numerical errors 
in the coefficient. Thus, we calculate the coefficients of H by the so-called least squares 
division [TH], followed by correcting the coefficients in F and G by using the calculated H, 
as follows. 

For polynomials F, G, A and B represented as in ([2]) and H represented as H{x) = 
{hd^i + hd^2i)x'^ + • ■ • + (/iQ.i + hf)^2i)x'^, solve the equations HB = F and HA = G with 
respect to H as solving the least squares problems of linear systems 

CdiA)*{hd,i + hd^2i,---,ho,i + /io,2i) = *(ffn,i + .9n,2i, • ■ • ,ff04 +.90,2*), (16) 

Gd{B)*{hd,i + hd^2i, ■ • ■ , ho,i + /io,2i) = *(/m,i + f,n,2i, . . . , /o,i + /o,2i), (17) 

respectively. Then, we transfer the linear systems (TTBl) and (flT)) . as follows. For (fT7|) . let us 
express the matrices and vectors as the sum of the real and the imaginary part of which, 
respectively, as Gd{B) = Bi + iB2,*(/id,i + hd,2i, ■ ■ ■ ,hQA + ^0,2*) = hi + i/i2,*(/m,i + 
/m,2*, . . . , /o,i + /o,2i) = /i + if2- Then, ^ is expressed as 

{Bi + iB2){hi + ih2) = {fi + if2). (18) 

By equating the real and the imaginary parts in Eq. (|18p . respectively, we have {Bihi — 
B2h2) = fi, {Bih2 + B2hi) = /2, or 



Bi -B2\ fhi\ ^ (h 
B2 Bi I \h2 \h 



(19) 



Thus, we can calculate the coefficients of H{x) by solving the real least squares problem 
([T5)) . We can solve (IT51) similarly. 

Let Hi{x), H2{x) G C[a;] be the candidates for the GCD whose coefficients are calculated 
as the least squares solutions of (fTO]) and (flT)) . respectively. Then, for i = 1,2, calculate the 
norms of the residues as r^ — \\F ~ HiB\W + \\G — HiAW"^, respectively, and set the GCD 
H{x) be Hi{x) giving the minimum value of r^. 

Finally, for the chosen H{x), correct the coefficients oi F{x) and G{x) as F{x) — H{x) ■ 
B{x), G{x) — H{x) ■ A{x), respectively. 



4 Experiments 

We have implemented the GPGCD algorithm for polynomials with the complex coefhcients 
on the computer algebra system Maple and compared its performance with a method based 
on the structured total least norm (STLN) method [7] for randomly generated polynomials 
with approximate GCD. The tests have been carried out on Intel Core2 Duo Mobile Pro- 
cessor T7400 (in Apple MacBook "Mid-2007" model) at 2.16 GHz with RAM 2GB, under 
MacOS X 10.5. 

In the tests, we have generated random polynomials with GCD then added noise, as 
follows. First, we have generated a pair of monic polynomials Fq{x) and Go{x) of degrees 
m and n, respectively, with the GCD of degree d. The GCD and the prime parts of 
degrees m — d and n ~ d are generated as monic polynomials and with random coefficients 
c G [—10, 10] of floating-point numbers. For noise, we have generated a pair of polynomials 
F^{x) and G-^{x) of degrees m — \ and n — 1, respectively, with random coefficients as the 
same as for Fq{x) and Gq{x). Then, we have defined a pair of test polynomials F{x) and 
G{x) as 

F{x) = Fo{x) + J Fn(x), G{x) = Go{x) + — ^£— Gn(x), 

respectively, scaling the noise such that the 2-norm of the noise for F and G is equal to ep 
and cg, respectively. In the present test, we set ep — ec = 0.1. 

In this test, we have compared our implementation against a method based on the 
structured total least norm (STLN) method [7], using their implementation (see Acknowl- 
edgments). In their STLN-based method, we have used the procedure C_con_mulpoly 
which calculates the approximate GCD of several polynomials in C[x]. The tests have been 
carried out on Maple 12 with Digits=15 executing hardware floating-point arithmetic. For 
every example, we have generated 100 random test polynomials as in the above. In execut- 
ing a modified Newton method, we set a threshold of the 2-norm of the "update" vector in 
each iteration e = 1.0 x 10~*; in C_con_mulpoly, we set the tolerance e = 1.0 x 10~*. 

Table [T] shows the results of the test: m and n denotes the degree of a pair F and G, 
respectively, and d denotes the degree of approximate GCD. The columns with "STLN" 
are the data for the STLN-based method, while those with "GPGCD" are the data for the 
GPGCD method. "Error" is the perturbation \\F - F\\l + \\G-G\\l, where "ae-6" denotes 
a X 10^**; "^Iterations" is the number of iterations; "Time" is computing time in seconds. 

We see that, in the most of tests, both methods calculate approximate GCD with almost 
the same amount of perturbations, while the GPGCD method runs much faster than STLN- 
based method by approximately from 10 to 30 times. Note that, in contrast to the real 
coefficient case [TTl, both methods have converged in all the test cases with the number of 
iterations and sufficiently small amount of perturbations as approximately equal to those 
shown as in Table [TJ 

5 Concluding Remarks 

Based on our previous research [17], we have extended our GPGCD method for polynomials 
with the complex coefficients. 

Our experiments have shown that, as in the real coefficients case |17j . our algorithm 
calculates approximate GCD with perturbations as small as those calculated by methods 
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Time (sec.) 
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STLN 
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STLN 


GPGCD 


1 


10,10 


5 


3.72e-3 


3.72e-3 


4.48 


4.43 


1.79 


0.15 


2 


20,20 


10 


4.16e-3 


4.16e-3 


4.24 


4.22 


5.88 


0.30 


3 


30,30 


15 


4.33e-3 


4.33e-3 


4.54 


4.48 


14.29 


0.58 


4 


40,40 


20 


4.48e-3 


4.48e-3 


4.08 


4.08 


24.10 


0.88 


5 


50,50 


25 


4.63e-3 


4.64e-3 


4.05 


4.12 


39.19 


1.36 


6 


60,60 


30 


4.61e-3 


4.61e-3 


4.02 


4.06 


60.48 


1.96 


7 


70,70 


35 


4.82e-3 


4.82e-3 


3.90 


4.02 


84.51 


2.66 


8 


80,80 


40 


4.84e-3 


4.84e-3 


3.88 


4.04 


116.03 


3.65 


9 


90,90 


45 


4.79e-3 


4.79e-3 


3.85 


4.01 


151.27 


4.66 


10 


100,100 


50 


4.77e-3 


4.78e-3 


3.83 


4.06 


199.48 


6.00 



Table 1: Test results for large sets of polynomials with approximate GCD. See Section H] 
for details. 

based on the structured total least norm (STLN) method, while our method has shown sig- 
nificantly better performance over the STLN-based methods in its speed, by approximately 
up to 30 times, which seems to be sufficiently practical for inputs of low or moderate de- 
grees. This result shows that, in contrast to their structure preserving method, our simple 
method can achieve accurate and efficient computation as or more than theirs in calculating 
approximate GCDs. 

Our future research includes theoretical investigation of convergence properties, inves- 
tigation for efficient computation in solving a linear system in each iteration by analysis of 
the structure of matrices, generalization of our method to several input polynomials, and 
so on. 
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